The main data used in this tutorial and in the lecture are about the geolocalisation of french restaurants in Paris and in a department called Haute-Garonne. We use two different sources:
SIRENE has the advantages of being rigorous and exhaustive on the French territory.
OSM has many benefits, ensuring transparent data provenance and ownership, enabling real-time evolution of the database and, by allowing anyone to contribute, encouraging democratic decision making and citizen science.
sf::st_read function also work? Why?
readRDS function.
st_read would not work because ‘iris_31.rds’ is not a shapefile but a file already R formatted. Simply load it with the readRDS function.
plot(iris_31). What do you notice ?
sf::st_geometry function? What solution do you propose then?
sf::st_geometry makes it possible to isolate the information contained in the ‘geometry’ column of the sf object. Using it, we put aside other variables (here CODE_IRIS, P14_POP and INSEE_COM).
crs="+proj=aeqd +lat_0=90 +lon_0=0” to see a clear difference and create a layer called ‘iris_31_aeqd’.
sf::st_crs function to guess the projection and sf::st_transform to change the projection.
Units: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00 51153.341 56509.22 53137.481 13756.98
[2,] 51153.34 0.000 20957.21 3086.011 49828.91
[3,] 56509.22 20957.212 0.00 11077.630 61863.18
[4,] 53137.48 3086.011 11077.63 0.000 54345.46
[5,] 13756.98 49828.910 61863.18 54345.458 0.00
Units: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00 57206.922 62306.32 59385.205 13798.28
[2,] 57206.92 0.000 20957.48 3204.855 55368.23
[3,] 62306.32 20957.481 0.00 11076.220 66744.65
[4,] 59385.20 3204.855 11076.22 0.000 59843.30
[5,] 13798.28 55368.225 66744.65 59843.302 0.00
[1] FALSE
dplyr package: select, group_by et summarize. These functions also work with sf objects.
cartography package, simply plot a map of french intercommunality with a proportional circle layer related to the number of restaurants.
propSymbolsLayer function allows you to draw proportional circles.
We would like here to design EPCI maps that combine the number of restaurants and the number of restaurants per 10,000 inhabitants.
Data preparation:
getBreaks et carto.pal functions of the cartography package. For the creation of the typo variable, you can use the cut function and apply the parameters digit.lab = 2 and include.lowest = TRUE.
library(sf)
library(cartography)
library(dplyr)
# Import data
fra <- st_read("../data/fra.shp", quiet = TRUE)
epci_31 <- readRDS("../data/epci_31.rds")
# Create the variable
epci_31$nb_rest_10000inhab <- 10000 * epci_31$nb_of_rest / epci_31$P14_POP
# Define breaks
bks <- getBreaks(v = epci_31$nb_rest_10000inhab, method = "quantile", nclass = 4)
# Define color palette
cols <- carto.pal("orange.pal", length(bks)-1)
# Create a "typo"" variable
epci_31 <- epci_31 %>%
mutate(typo = cut(nb_rest_10000inhab,breaks = bks, dig.lab = 2,
include.lowest = TRUE))cartography package, make the following map which contains in a choropleth layer the variable nb_rest_10000inhab and in a proportional circle layer the variable nb_of_rest. Do the same thing with the ggplot2 package.
With cartography:
# Define plot margins
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
# Find EPCI bounding box
bb <- st_bbox(epci_31)
# Plot France using EPCI boundingbox
plot(st_geometry(fra), col="ivory", border = "ivory3",
xlim = bb[c(1, 3)], ylim = bb[c(2, 4)])
# Plot the choropleth layer
choroLayer(epci_31, var = "nb_rest_10000inhab",
breaks = bks, col = cols, border = "grey80", lwd = 0.5,
legend.pos = "topleft",add = TRUE,
legend.title.txt = "Number of restaurants\nfor 10,000 inhabitants")
# Plot proportionnal symbols
propSymbolsLayer(epci_31, var="nb_of_rest", col="#440170",border=NA,
legend.pos="left", inches=0.4, add = TRUE,
legend.title.txt = "Number of restaurants")
# Add a layout layer
layoutLayer(title = "Restaurants", sources = "Insee, 2018",
author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred",
coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
# Add a north (south) arrow
north(pos = "topright", south = TRUE)With ggplot2:
library(ggplot2)
map_ggplot <- ggplot() +
geom_sf(data = fra, colour = "ivory3",
fill = "ivory") +
geom_sf(data = epci_31, aes(fill = typo), colour = "grey80") +
scale_fill_manual(name = "Number of restaurants\nfor 10,000 inhabitants",
values = cols, na.value = "#303030")+
geom_sf(data = epci_31 %>% st_centroid(),
aes(size= nb_of_rest), color = "#440154CC", show.legend = 'point')+
scale_size(name = "Number of restaurants",
breaks = c(1, 500, 3200),
range = c(0,18))+
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(epci_31)[c(1,3)],
ylim = st_bbox(epci_31)[c(2,4)]
) +
theme_minimal() +
theme(panel.background = element_rect(fill = "azure",color=NA)) +
labs(
title = "Restaurants",
caption = "Insee, 2018\nKim & Tim, 2018"
)
plot(map_ggplot)cartography package.
propSymbolsChoroLayer function allows you to draw colored proportional circles.
# Define plot margins
par(mar = c(0.2, 0.2, 1.4, 0.2), bg = "azure")
# Find EPCI bounding box
bb <- st_bbox(epci_31)
# Plot France using EPCI boundingbox
plot(st_geometry(fra), col="ivory", border = "ivory3",
xlim = bb[c(1, 3)], ylim = bb[c(2, 4)])
# Plot EPCI
plot(st_geometry(epci_31), col="ivory3", border = "ivory2", add=T)
# Plot the choropleth layer
propSymbolsChoroLayer(epci_31, var = "nb_of_rest", var2 = "nb_rest_10000inhab",
breaks = bks, col = cols, border = "grey80", lwd = 0.5,
legend.var.pos = "topleft", legend.var2.pos = "left",
add = TRUE, inches = 0.4,
legend.var.title.txt = "Number of restaurants",
legend.var2.title.txt = "Number of restaurants\nfor 10,000 inhabitants")
# Add a layout layer
layoutLayer(title = "Restaurants", sources = "Insee, 2018",
author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred",
coltitle = "white", postitle = "center",
frame = TRUE, scale = 10)
# Add a north (south) arrow
north(pos = "topright", south = TRUE)cartography package, display the number of restaurants and the number of restaurants per 10,000 inhabitants at the municipalities and EPCI scales. The two maps displayed side by side should be as much comparable as possible.
library(sf)
library(cartography)
library(dplyr)
# Import data
fra <- st_read("../data/fra.shp", quiet = TRUE)
epci_31 <- readRDS("../data/epci_31.rds")
com_31 <- readRDS("../data/com_31.rds")
# Create the variable
epci_31$nb_rest_10000inhab <- 10000 * epci_31$nb_of_rest / epci_31$P14_POP
com_31$nb_rest_10000inhab <- 10000 * com_31$nb_of_rest / com_31$P14_POP
# Define breaks for municipalities (we will use the same breaks for both maps)
bks_com <- getBreaks(v = com_31$nb_rest_10000inhab[com_31$nb_of_rest>0],
method = "quantile", nclass = 6)
# Define color palette
cols <- carto.pal("wine.pal", length(bks_com)-1)
# Define plot margins
par(mar = c(0, 0.1, 1.2, 0.1), bg = "azure", mfrow = c(1,2))
# Find EPCI bounding box
bb <- st_bbox(epci_31)
# Plot France using EPCI boundingbox
plot(st_geometry(fra), col="ivory", border = "ivory3",
xlim = bb[c(1, 3)], ylim = bb[c(2, 4)])
# Plot EPCI
plot(st_geometry(epci_31), col="ivory3", border = "ivory2", add=T)
# Plot the choropleth layer
propSymbolsChoroLayer(epci_31, var = "nb_of_rest", var2 = "nb_rest_10000inhab",
breaks = bks_com, col = cols, border = "ivory3",lwd = 0.6,
legend.var.pos = "bottomright", legend.var2.pos = "n",
add = TRUE, inches = 0.5,
legend.var.title.txt = "Number of restaurants")
# Add a layout layer
layoutLayer(title = "Restaurants", sources = "Insee, 2018",
author = "Kim & Tim, 2018",
theme = "green.pal", col = "darkred",
coltitle = "white", postitle = "center",
frame = FALSE, scale = NULL)
bb <- st_bbox(epci_31)
# Plot France using EPCI boundingbox
plot(st_geometry(fra), col="ivory", border = "ivory3",
xlim = bb[c(1, 3)], ylim = bb[c(2, 4)])
# Plot EPCI
plot(st_geometry(com_31), col="ivory3", border = "ivory2",lwd = .5, add=T)
# Plot the choropleth layer
propSymbolsChoroLayer(com_31, var = "nb_of_rest", var2 = "nb_rest_10000inhab",
breaks = bks_com, col = cols, border = "ivory3",lwd = 0.6,
fixmax = max(epci_31$nb_of_rest),
legend.var.pos = "n", legend.var2.pos = "bottom",
add = TRUE, inches = 0.5, legend.var2.values.rnd = 0,
legend.var2.title.txt = "Number of restaurants\nfor 10,000 inhabitants")
# Add a layout layer
layoutLayer(title = "Restaurants", sources = "",
author = "",
theme = "green.pal", col = "darkred",
coltitle = "white", postitle = "center",
frame = FALSE, scale = 10)
# Add a north (south) arrow
north(pos = "topright", south = TRUE)mapview package. Try using different parameters to customize your map.
map.types, col.regions, label, color, legend, layer.name, homebutton, lwd … parameters of the mapview function.
library(mapview)
library(sf)
library(cartography)
sir_31 <- readRDS("../data/sir_31.rds")
mapview(sir_31, map.types = "OpenStreetMap",
col.regions = "#940000",
label = paste(sir_31$L2_NORMALISEE, sir_31$NOMEN_LONG, sep = " - "),
color = "white", legend = TRUE, layer.name = "Restaurants in SIRENE",
homebutton = FALSE, lwd = 0.5)